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ABSTRACT 

Numerical simulations of galaxy clusters including two species - baryonic 
gas and dark matter particles - are presented. Cold Dark Matter spectrum, 
Gaussian statistics and flat universe are assumed. The dark matter component 
is evolved numerically by means of a standard particle mesh method. The 
evolution of the baryonic component has been studied numerically by using a 
multidimensional (3D) hydrodynamical code based on modern high resolution 
shock capturing techniques. These techniques are specially designed for treating 
accurately complex flows in which shocks appear and interact. With this 
picture, the role of shock waves in the formation and evolution of rich galaxy 
clusters is analyzed. Our results display two well differenced morphologies of 
the shocked baryonic matter: filamentary at early epochs and quasi- spherical at 
low redshifts. 

Subject headings: cosmology: theory — hydrodynamics — large-scale structure of 
the universe — methods:numerical — shock waves 
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1. Introduction 

Galaxy clusters are the largest systems gravitationally bounded in the Universe. Their 
study has been a fashion topic in Cosmology since last years. Work on topics related 
with galaxy clusters is worthly to: i) understand the formation, evolution, dynamics and 
morphology of these systems, ii) learn on the physical processes involved in them, and 
iii) find out some information concerning with fundamental parameters in Cosmology as 
density parameter , Bubble's constant (if), and the spectrum of the primordial density 
field. 

During last years technical improvements have produced huge quantity of data about 
galaxy clusters. Let us mention the new galaxy surveys (Guzzo 1996, and references 
therein) and the extensive observation in X-rays using sateUites as ROSAT or ASCA. These 
huge volume of data strongly motivates a lot of theoretical work trying to explain the 
observational results. 

Prom the theoretical point of view, numerical simulations are the best tools to 
understand physics involved in galaxy clusters. At the beginning, numerical simulations 
of galaxy clusters where performed using N-body techniques. Since then, they have been 
extensively used and have produced important results ( see, e.g, Efstathiou et al. 1985, 
Bertschinger & Gelb 1991, Xu 1995). Next step in the full description of galaxy clusters was 
to introduce in the picture a baryonic component. The numerical methods developed in 
order to deal with baryonic matter were more sofisticated and expensive in computational 
resources. As a consequence, it was not possible to carry out numerical simulations with 
two species (dark matter and baryonic gas) until late eighties. 

Cosmological hydrodynamic codes have been usually classified in two main categories: 
a) the so-called Lagrangian methods, like the Smoothed Particles Hydrodynamics (SPH) or 
ulterior extensions based on them, and b) Eulerian codes. 
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SPH methods were first proposed by Gingold & Monaghan (1977), and Lucy (1977). 
Among the best features of this technique, it should be pointed out its high resolution in 
dense regions. This property is directly derived from its Lagrangian character. The first 
implementations of SPH techniques had some weak points: i) The low density regions were 
badly described due to the Lagrangian character of the method, ii) Discontinuities and 
strong gradients were poorly solved and an important diffusion was introduced, iii) They 
were not conservative. Nevertheless, these previous problems were overcome in the modern 
implementations of these techniques. Improved SPH techniques have been widely developed 
for cosmological applications (see, e.g., Evrard 1988, Hernquist & Katz 1989, Navarro & 
White 1993, Gnedin 1995). 

Numerical cosmological codes using an Eulerian approach to study baryonic gas inside 
galaxy clusters have been also developed. Some of these hydro-codes use artificial viscosity 
in order to deal with shock waves (Cen 1992, Anninos et al. 1994). These techniques require 
a good calibration of the free parameters which are introduced by hand and state some 
numerical problems. Recently, a new family of finite difference methods, which use Eulerian 
approaches and avoid artificial viscosity, has been developed in numerical Cosmology. They 
are the so-called high resolution shock capturing methods (HRSC), the modern extensions 
of the original Godunov's idea (1959). According to the Riemann solver and the procedure 
in order to achieve spatial accuracy, we can distinguish three groups: 1) the ones following 
Harten's scheme (1983), like Ryu et al. (1993), 2) those using the analytical solution of 
the Riemann problem for the Newtonian dynamics of ideal gases and the PPM scheme 
described by CoUela & Woodward (1984) , like Bryan et al. (1994), and 3) the codes using 
Roe's Riemann solver (Roe 1981) plus the MUSCL or PPM cell reconstruction, like in 
Quilis et al. (1996). In this last reference, the code used in present paper is described and 
tested appropriately. 
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An exhaustive comparison among all these kinds of cosmological hydrodynamic codes 
can be found in Kang et al. (1994). 

Due to the Eulerian character of our code, it does not show -in dense regions- a 
resolution as good as the Lagrangian ones, and it requires more computational resources. 
However, HRSC schemes -by construction- have excellent properties in order to deal with 
shocks, discontinuities, and strong gradients. HRSC techniques typically solve shocks in 
two cells. Due to their intrinsic properties, the detection of shocks is independent on the 
number of cells used in the simulations. It should be pointed out that this last property is 
really important when three-dimensional simulations are carried out. In these simulations 
the size of the grid is a stringent constraint due to its high cost in computational resources. 
Moreover, these methods are conservative by construction, that is , quantities which should 
be physically conserved are numerically conserved up to the order of the method. It should 
be also noticed that these methods show good results in extreme low density regions 
(Einfeldt et al. 1991). 

As it has been pointed out by several authors, the role of shock waves can be extremely 
important in order to understand the heating processes in the intracluster medium (ICM). 
In this paper wc are interested in understanding and quantifying the role of shocks. In 
order to do that it is crucial to use numerical codes able to manage with complex flows. 
Yes, indeed, one of the important features of HRSC techniques is just to treat numerically 
shocks and strong discontinuities giving sharp profiles (in a few numerical cells, as we have 
mentioned above) independently of the size of the grid. Hence, formation, evolution, and 
interaction of shocks in 3D flows can be analyzed accurately with HRSC schemes, and, 
consequently, their use is absolutely justified in order to study shocks and their consequences 
on the ICM's dynamics. 

Hereafter, t stands for the cosmological time, to is the age of the Universe, a{t) is the 



- 6 - 



scale factor of a flat background. Function d/a is denoted by where the dot stands for 
the derivative of a with respect to the cosmological time. Hubble constant is the present 
value of if; its value in units of 100 Km Mpc~^ is the reduced Hubble constant h. In 
our computations we have assumed h = 0.5 . Velocities are given in units of the speed of 
light. Baryonic, dark matter, and background mass density are denoted by p^, p^^^^, and p^, 
respectively. The density contrast is 6^ = {pb ~ Pb)/Pb ^^r baryonic matter, analogously is 
defined 5^^^,^ for dark matter. 

The plan of this paper is as follows: In Section 2, our numerical cluster model is 
described. In Section 3, the results of the simulations are analyzed. Finally, a general 
discussion is presented in Section 4. 



2. Cluster Model 

2.1. Initial conditions 

Initial conditions are given at redshift z = 100. We have considered a fiat universe 
{Qo = 1) and a CDM scenario. The density profile is obtained by using the power spectrum: 

\h\' = (1) 

{1 + (3k + ukl + -fk^y 

where p = 1.7 (noh^)-^ Mpc, u = 9.0{noh^)~^-^ Mpc^-^, 7 = {VLoh'^)-'^ Mpc^ , and the 
constant A (normalization) has been fitted by using erg = 0.63, where ag = ((^)^)^ on the 
8 h^^ Mpc scale. Assuming that the density field is Gaussian, realizations in the position 
space can be obtained from Eq.(|l|). 

Hoffman and Ribak (1991) described a method to produce realizations of Gaussian 
fields under some constraints. This procedure was extended by van den Weygaert and 
Bertschinger (1996). We have used this powerful method in order to obtain the initial 
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conditions for our cluster. In particular, we have assumed three constraints: 1) the cluster 
peak is located just at the center of the box, its amplitude corresponds to the maximum 
density contrast in our initial box (64 Mpc per edge in comoving coordinates), 2) the mass 
enclosed inside a Gaussian ball, M — {27:^^ ^'^^pbR} with Rf — Ah~^ comoving Mpc, is 
M = 6 X 10^^ Mq -centered at the peak position- , and 3) the amphtude of the initial 
density contrast -cluster peak- is 3cr after smoothing on scales ~ ih'^Mpc using a Gaussian 
window function. 

The 3D Zeldovich's approximation is used to evolve from the above initial conditions 
up to 2; = 30. At this redshift, a comoving box of 20 Mpc per edge is extracted (zoom) from 
the initial box. The maximum of the density contrast is located at the center of the new 
box. 

The density and the velocity profiles of the initial perturbation are known after 

performing the zoom ai z = 30. Then, we define the composition of this perturbation, 

namely, we fix the ratio between the dark matter (DM) and baryonic densities. It is 

considered that the 90 % of the matter is dark and the remaining is baryonic. The DM 

component is discretized into particles of the same mass (the mass of the particles depend 

upon the resolution of each simulation). The baryonic component is considered as an ideal 

fluid of adiabatic exponent 7, fully ionized, with mass abundances of 75% Hydrogen and 

25% Helium, which represents the intracluster medium (ICM). We assume that the initial 

velocities of both components are identical. The initial homogeneous temperature is set up 

assuming that at 2; = 30 the temperature is given by T30 = T200 \ -^—7 ?r2r 1 ) where 

\Pb{z^ 200) J 

T200 — 2.73(1 + z) with z — 200, and 7 is the adiabatic exponent (see , e.g., Anninos & 
Norman 1996). 

The above initial conditions give the seed for the formation of a large Abell cluster. 
Some comments about the boundary conditions are deserved. The baryonic and DM 
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have low densities at the faces of the numerical box. Moreover, DM particles are not allowed 
to escape out the box. In order to do that, when one particle crosses some of the faces of 
the computational box, it is introduced through the symmetric face with the same velocity. 
Thus, a reasonable approximation to periodic boundary conditions is considered. 



2.2. Nonlinecir evolution 

The nonlinear evolution of the baryonic component is described by the hydrodynamics 
equations (Peebles 1980): 

^ + iv.(l + .>=0 (2) 
dv 1 11 

^ + -(v ■ V)v + Hv^ Vp--V0 (3) 

at a Pf^a a 

^ + - V • [(£; + p)v\ = -m{E + p) - Hp.if' - ^ V0 - A (4) 



and the evolution of the DM particles obeys the particle equations: 

dx V 
dt a 

dv _ V0 
dt a 



(5) 

- Hv (6) 



where x, v — = {'t^x,'t^y,t^z), E and <f){t,x) are, respectively, the Eulerian coordinates, 

the peculiar velocity, the total energy E — p^€+ \p^v'^ {v"^ = '^^ + '^^ + '^D ' pecuhar 
Newtonian gravitational potential. The pressure and the internal energy per unit mass are 
denoted, respectively, as p and e. The above system of hydrodynamical equations is closed 
with the equation of state of an ideal gas p— (7 — l)p^e and 7 = 5/3 (monoatomic gas). 

In some applications described in the present paper two cooling processes are taken 
into account, Compton coohng (A'-^) and thermal Bremsstrahlung (A'^''). The well known 
expressions for these processes are (see, for instance, Umemura & Ikeuchi (1984)): 

A^ = 5.4 X 10-^^(1 + zYueT (ergcm-^s-^) (7) 
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and 

A^'^ = 1.8 X IQ-^'^neT^nHH + ^UHeiii) {ergcm-^s-^) (8) 

where ne,nHn, and UHeiii are the number density of electrons, protons, and hehum nuclei, 
respectively. The cooling effects are introduced in the source term A = A*^ + A^** in Eq.(^). 

Both components, gas and particles, are gravitationally coupled through Poisson's 
equation: 

= lH'a'6T, (9) 

where 6t = {p,, + P^m ~ Pb)/Pb total density contrast. 

The hydrodynamical equations define a hyperbolic system of conservations laws 
(with sources). This important property allows us to apply a new powerful family of 
numerical methods (HRSC) in order to solve this system. The main features of these 
methods are summarized: 1) the quantities which should be physcially conserved, are 
numerically conserved by construction, due to the fact that the algorithm is written in 
conservative form , 2) these methods are of high order in the smooth regions of the flow, 
3) shocks are sharply solved in typically two numerical cells, independently of the number 
of cells used in the simulations , 4) numerical artifacts like artificial viscosity are avoided, 
5) these methods are free of spurious oscillations around shocks, and 6) they are able to 
handle strong gradients and discontinuities. In Quilis et al. (1996) we took advantage of the 
hyperbolicity property of this system in order to design a multidimensional hydrodynamic 
code. A full description of this code and a set of tests proving its performance were 
presented in this reference. The above method is used in present paper in order to study 
the evolution of the baryonic gas. 

The evolution of DM particles described by Eqs. (^-|§) is studied through a standard 
Particle Mesh (PM) method (Hockney & Eastwood 1988). This method has some important 
features from the point of view of our applications: 1) the PM methods are widely extended 
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and well known, 2) they are easily programmable and can be directly coupled with a 
hydro-code, and 3) their CPU cost is quite low in comparison with others particle type 
methods. This last property is crucial when the particle code must be coupled with a 
hydro-code like the one described above where the computational effort is severe. 

Poisson's equation is solved by using Fast Fourier Transform (FFT) methods 
(Press et al. 1987). In order to recover a continuous density contrast field for dark matter 
component, we use a standard cell-in-cloud (CIC) scheme (Hockney & Eastwood 1988) at 
each time step. 

The time step is computed following Quilis et al. (1996). We have considered Courant's 
and dynamical characteristic times. These times must be corrected in order to ensure the 
numerical stability of the algorithm. Typical values of these correction factors, known 
as CFL, are CFLcourant = 0.6 and CFLdynamicai ~ 10~^. Besides the characteristic times 
mentioned before, cooling processes introduce a new characteristic time. Nevertheless as 
it is known (Tsai et al. 1994), for the scales considered in this paper, the characteristic 
cooling time is larger than the dynamical one and, consequently, cooling would not expect 
to be an effective process. 

In order to display the conservation properties of the whole code - hydro, particles 
and the Poisson solver - the global energy conservation is investigated in the Appendix A. 
Some comments on the numerical spatial resolution of our simulations are also discussed in 
Appendix B. 

3. Results 

In this Section we are going to discuss the results obtained running our code in two 
cases: i) without cooling processes (denoted by A), and ii) including cooling effects (denoted 
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byB). 

Numerical simulations have been carried out using 128^ DM particles and 128^ 
baryonic cells. The physical comoving box size is 20 Mpc per edge. Since we are interested 
in the study of ICM and not concerned with smaller scales like galaxies theirselves, the 
choice of 128^ DM particles and 128^ baryonic cells seems adequate. The resolution of the 
hydrodynamical code is ~ 0.15 comoving Mpc, however the PM code resolutions is two 
cells. Hence, the minimal resolution of our simulations will be around ~ 0.3 comoving Mpc. 
Taking into account the properties of HRSC techniques, the grid used in present paper 
seems a reasonable compromise between the resolution needed in order to describe physical 
processes and our current computational limitations. 

The CPU cost of our numerical simulations is 183.2 s(25.1 s) per time step (with 128^ 
numerical cells, and 128^ particles) in a Hewlett-Packard J280 (Silicon Graphics Origin 
2000 with 8 processors for parallel version). 

3.1. Cluster evolution 

In the standard CDM scenario, clusters of galaxies grow by merging from smaller 
collapsed structures. This behaviour can be tested through N-body simulations . Figure 1 
shows the DM particle distribution for our A simulation at six epochs. DM dynamics seems 
to be not much affected by gas dynamics and cooling effects. By this reason the results 
from B simulations are not shown. In Figure 1, 32^ particles have been plotted. The edges 
of the boxes are 20 comoving Mpc length. 

The implementation of a baryonic gas complicates the description of the cluster. 
This is a component with pressure and with a different dynamics. This new component 
could undergo compression, expansion, heating and cooling. Non-adiabatic processes 
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like shock waves and the coohng - due to Bremsstrahlung and Compton - could also be 
present. These deviations from adiabaticity are very important in order to understand the 
mechanism which provides an extra energy to ICM. Figure 2 shows ICM density contrast 
colour contours for A simulation. The slices are 20 (0.15) comoving Mpc length (depth). 
Four times are displayed. The evolutionary picture starts with collapsing processes of 
small substructures. These substructures undergo a merging process towards a well defined 
central object. This picture is valid for both simulations, A and B. Since no significant 
differences can be found, for the sake of briefness only A results have been displayed. These 
results supports the well known fact that for the scales considered in this paper, the cooling 
is not a relevant process (see ,e.g, Tsai et al. 1994). 

A good description of the cluster evolution is given by the profiles, at different epochs, 
of the main quantities averaged in spherical shells. These shells are located at several radii 
from the cluster center. This center is defined as the minimum of the total gravitational 
potential. Each shell has a fixed comoving width. An average value inside the shell can 
be estimated for each variable. Figure 3 plots the average density contrast for baryonic 
gas (top panel) and DM (bottom panel) versus radial distance to the cluster center in 
comoving Mpc. In each panel several curves are shown corresponding to different times. 
As no significant differences arise between A and B simulations, only case A is presented. 
It should be pointed out the resemblance of the profiles plotted here and the ones shown 
by Navarro et al. (1995). As it was explained in detail in this reference, for the spatial 
resolution considered in this paper, both components evolve keeping the original ratios 
between their densities. 

The gross properties of the cluster dynamics can be fairly inferred studying average 
radial velocities of DM and baryonic gas. These averages are computed in a similar way to 
the density contrast. Figure 4 plots the average radial velocities, in units of c, for baryonic 
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gas (top panel) and DM particles (bottom panel). As in Figure 3, only A results are shown. 

From Figs. 3 and 4 the following stages in the evolution could be roughly distinguished. 
First, an inf ailing phase until 2; ~ 2, characterized by a typical profile in velocity (gas 
and DM) and increasing in density contrast (both components). Let us point out that 
the profile, during the infalling epoch, of the baryonic component remind us the typical 
one in any collapsing object, the minimum of which separates a subsonic homologously 
part (the inner one) from the supersonic freely falling part (the outer one). In a second 
phase, gas velocities change their signs in the inner part while outer shells go on falling. 
Gas density contrast slightly decreases due to the gas outfiow. Nevertheless, this effect is 
subdominant, and the ratios between the averaged density profiles of gas and DM undergo 
small variations. DM particles are not affected by this dynamics and exhibit some tendency 
to reach a dynamical equilibrium. 

3.2. Looking for shocks 

In the present Section we are specially interested in giving a set of criteria allowing us 
to identify those regions where shocks form and evolve. In the next Section will analyze the 
energy released to the ICM by these shocks. 

From the top panel of Figure 4 some strong evidences of the existence of a quasi- 
spherical shock are pointed out. The velocity profiles show how the evolution begins with a 
typical infalling profile. This velocity is strongly correlated with the growing of the baryonic 
density contrast, for the same time, in top panel of Figure 3. As the density rises , the 
pressure becomes more important. At 2; ~ 2, the velocities at the inner part of the cluster 
change their sign, and the gas moves out. This is an evidence of a quasi-spherical shock 
moving outwards from the core cluster. 
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FoUowing Anninos and Norman (1996) we study the entropy s — ln{T/p2~^). In 
an adiabatic process this quantity should keep constant. Variations in s are evidence of 
non-adiabatic processes. In particular shock waves strongly influence the variations in the 
entropy. The maximum in the entropy could identify the position of shocks. Figure 5 plots 
the average entropy in spherical shells for A case at several times. 

The average entropy profiles evidence the existence of a quasi-spherical shock travelling 
outwards from the cluster center. The formation and evolution of the local maxima in 
the average entropy profiles can be directly correlated with the changes in the velocity 
profiles for baryonic gas (top panel Figure 4). The averaging procedure in shells is itself 
an important source of numerical diffusion. When the shock deviates from sphericity or 
there are a set of smaller and irregular shocks, then the average in spherical shells does 
not produce two regions separated by a sharp jump in the entropy profile. In Figure 5 
two regions are clearly differenced. The pre-shock region has a roughly constant entropy, 
while the post-shock shows a decreasing inward entropy with values higher than those of 
pre-shock region. 

The above analysis evidences the existence of, at least, one quasi-spherical shock. 
Nevertheless, it is not possible to know if there is a unique quasi-spherical shock, or there 
are some shocks at several scales which in average give the global shock in Figure 5. From 
this Figure, no ideas about the scales and form of the shocks can be obtained. In order 
to clarify this point. Figure 6 shows slices for the divergence of baryonic velocity, entropy, 
temperature and pressure. These slices are centered at the point where the gas density is 
maximum, and are projected along z-axis. Size and depth of the slices are the usual through 
this paper (see Fig. 2). 

As it has been mentioned above, entropy, pressure and temperature exhibit a well 
known behaviour in shocks. Nevertheless, some comments on the usefulness of studying 
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the divergence of baryonic velocity are worthwhile. If we assume that, at a given time, the 
density contrast (5f, is homogeneous, and that the source terms in Eq. are negligible, 
then, the time variation of the density contrast is related with the divergence of the velocity 
through the following equation, 

dln5. 1„ ^ 

— -A = — V-^ 10 
at a 

This equation points out two different regimes: 1) compression, i.e., V ■ w < 0, and 2) 
rarefaction, i.e., V ■ "U > 0. Therefore, the changes in the sign of the divergence of the 
velocity allow us to distinguish between compressed and rarefied regions. This discussion 
can always be applied locally, and it gives a roughly description of the dynamical state of 
the system. 

In our simulations, neighbouring regions might be separated by a shock if simultaneously, 
exhibits large changes , the entropy increases clearly, and the temperature and pressure 

show strong gradients. The features of the slices for these quantities, allow us to have a 

good description of the shock scales and their morphology. 

Strong evidences of a quasi-spherical shock travelling outwards appear in all the 
quantities plotted in Figure 6. A thin shell, where the divergence of baryonic velocity 
and temperatures change strongly and the entropy reaches their maxima values, is clearly 
appreciable. Pressure exhibits a tiny jump in the same region, but, unfortunately, due to 
the extreme variations in pressure from the center to the outer regions, the colour scale is 
not able to reproduce this pressure jump. 

Even more, it is possible to define some criteria which could label a cell which is 
involved in a shock. Following the criteria used by Colella & Woodward (1984) to detect 
shocks, we generalize them for a multidimensional case. A numerical cell, (i, j, /c), is labelled 
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as involved within a shock if the following conditions are satisfied: 



\Pi+l,j+m,k+n Pi—l,j—m,k—n\ 
Pi,j,k 

/ \ 9 1 



> Pi l,m,n = 0,1 (11) 



max 



W,j+^"^,k+n ^-l,j-m,k-n\ \ I , m, u = 0, 1 and q = X, y , z (12) 



\ f^^i-\-l .j-\-m,.k-\-rL P^i — l,j — m,k — n I 



> /33 Z,m,n = 0,1 (13) 



where /3i,/32 and /^s are three parameters to be fixed. Let us comment the above criteria to 
identify a cell as a shocked one. As it is well known in ID, across a shock, quantities such 
as p, f , and p have a jump discontinuity. The lack of numerical resolution may produce 
misleading results without a fine tuning of (3 parameters. For example, in the case of low 
values of cells where there are only strong gradients could be identified as shocked cells. 
These problems can be extended to the other conditions. Some special attention deserves 
strong rarefactions. Since regions where the gas is being fast rarefied could satisfy the above 
conditions even with high values for (3 parameters, one new condition must be considered 
in order to avoid identifying strong rarefactions as shocks. This condition takes advantage 
from the fact that the velocity in the post-shock region is always greater than that of the 
pre-shock region (Colella & Woodward 1984), i.e.. 



i~l,j,k > "^i+lj.fc (^"^ '^i,j-l,k > ^i,j+l,k ^'^ '^lj,k-l > ^lj,k+l i^^) 



V. 



Some experimentation has been carried out so as to get the values for the parameters 
/?i,/52, and /^s, and also to check their performance. As we have seen above, too low values 
for the parameters could be dangerous. However, high values could fail when the shocks 
present in the simulation were not too strong. In order to check conditions (|llHH) we need 



a simulation involving shocks, contact discontinuities and rarefactions, based on a problem 
with a well known analytical solution. The results of a 3D analytical test presented by 
Quilis et al. (1996) have been used. The identification of the cells involved in the shock is 
excellent in this case. Although general results are more complex, the method seems to be 
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able to identify shocked cells reasonably. In the practical implementation presented in this 
paper, /? parameters have been chosen under the constraints that the regions labelled in 
Fig. 6, as shock candidates, are most accurate recovered. The values used in this paper for 
those parameters, which fulfil the above criteria, are /3i = 3, /32 = 2, and (3^ — 2. 

Collecting together all the previous prescriptions. Figure 7 shows some 3D plots of the 
shocked cells at several redshifts for A simulation. At high redshifts, the shocked cells form 
some sort of filaments which would be following the ones arising in DM (see Fig. 1). Then, 
we could say that these shocks , which appears at quite early times, come from the collapse 
of the substructure and from merging processes among these substructures. Between 2; ~ 2 
and 2; ~ the shocked cells mainly trace a quasi-spherical shock moving outwards from 
the cluster center. Morphologies of the shocks display filamentary structures at the early 
epochs - high redshift - while at lower redshift a quasi-spherical shock appears at the inner 
regions and propagates outwards dominating the global structure. 

The presence of shocks is directly related with the dynamics of gas. By using the 
above conditions the number of shocked cells can be computed. Figure 8 plots the number 
of shocked cells versus time evolution for A(B) continuous (dashed) hne. The maximum 
number of cells involved in shocks is at 2; ~ 2. This epoch is directly correlated with those 
times when the quasi-spherical shock in Figs. 3 and 4 starts to form. 

3.3. Energy released by shocks 

In previous section the presence of shocks has been discussed. Now we try to quantify 
the energy released by those shocks to ICM. 

An easy way to quantify the shock heating is to compare some quantities obtained by 
the simulations with the values obtained considering an analytical adiabatic evolution. The 
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quantities derived from this analytical adiabatic evolution are labelled with tilde (~). As it 
is well known in the adiabatic regime the pressure and temperature are related with density 
at any time in the following form: 

P = i^pI (15) 

T = r^pr' (16) 

where k and t] are two constants. In order to fix the values of these constants we take the 
same values of pressure (p), density (p), and temperature (T) as the ones used, at the initial 
time, in the numerical simulations. 



If the simulated cluster would evolve pure adiabatically, Eqs. ([TB|-|TB|) should give. 



at any time, the same values of pressure and temperature than those obtained from the 
hydro-code. In particular, if the baryonic density derived from the code were used as an 
input for the above adiabatic laws, and the cluster would evolve pure adiabatically, no 
differences should be found, at any time, between results from Eqs. ( p!5| - |T6|) and the ones 
given by the numerical code. In practice, we compute all the pure adiabatic quantities 
{p, T, pb) using as input the density coming from the simulations. Any deviations between 
results from the simulation A, and results from Eqs. ([T5| - |T6D - where it has been used the 



density evolved by the A simulation - should be considered as a proof of the presence of 
non-adiabatic processes. 

The total internal energy for the whole box is defined as : 

Eu= [ PbedV (17) 
Jv 

being y the physical volume. 

Three total internal energies can be computed. First (second) is the one derived from 
the A(B) simulation. A third one is obtained from the analytical adiabatic evolution by 
using Eqs. (|K]-|T^). Top panel of Figure 9 shows total internal energy from A simulation 
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(continuous line) and B simulation (dashed-pointed line), and from analytical adiabatic 
evolution (dashed line). The curves from simulations appear clearly different from the 
analytical case. If the evolution were completely adiabatic, the internal energy for the A 
simulation (no cooling effects) should be identical to the one for the analytical adiabatic 
evolution. As this simulation does not consider cooling effects, all these differences are due 
to non-adiabatic processes. That means that shocks are the responsible of the most part 
of the ICM heating. The total internal energy for the B simulation is slightly lower than 
the one obtained in the A case at some phases of the evolution. These small differences are 
negligible in front of the differences between A or B cases against the analytical adiabatic 
evolution curve. Consequently, cooling effects are negligible compared with shocks. 

Total kinetic energy is estimated from the following volume integral: 



Bottom panel of Figure 9 plots time evolution of the ratio Eu/Ek for the data 
generated by our simulations. Continuous (dashed-pointed) line corresponds to A(B) 
simulation. These curves show those times of the evolution where the internal energy is 
being released to the ICM by shocks. The decreasing branch in Figure 8, after the first 
maximum, is related with the epoch in which shocks - already formed - are propagating 
and loosing energy, which contributes to the ICM heating. This branch is correlated with 
the raising branch of the Eu/Ek curve in Figure 9. 

Figure 10 is a plot of the time evolution of the logarithm of the average temperature in 
a ball of one comoving Mpc radius centered at the core cluster. As usual, data generated by 
A (continuous line) and B (dashed-pointed line) simulations, and the analytical adiabatic 
evolution (dashed line) are displayed. As the observations suggest, the values of temperature 
obtained in our simulations range around ~ 3 x 10^ K at low redshifts. Both simulations, 
A and B, show how the temperature has increased about six orders of magnitude over 
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that of the pure adiabatic evolution. There are no important differences between A and B 
temperatures, as in the case of the quantities analyzed previously. 

Figure 11 plots - for the case B - the X-Ray thermal bremsstrahlung luminosity 
produced by a ball of one comoving Mpc radius centered at the core cluster. The values are 
compatible with those given by the X-Ray observations which range between 10"^^ and 10^^ 
erg/s. 

4. Discussion 

In this paper we have used some numerical techniques recently applied to Cosmology. 
These techniques , HRSC, seems to be the most suitable in order to study the role of shocks 
in galaxy cluster evolution. The choice is justified by their properties to handle shocks. The 
capability of these techniques to capture shocks with very small diffusion is independent 
of the resolution used in the numerical simulations. Hence, by construction, shocks are 
captured even using coarse grids. This property is crucial in 3D applications. 

Previous sections illustrate the fact that non adiabatic processes, due to shocks, take 
an important role in the description of the ICM. In the model presented in this paper, that 
is, a baryonic fluid plus dark matter component coupled gravitationally, shocks are able to 
heat the ICM until values compatible with observational data. 

The calculations have been carried out in two cases: with and without cooling 
processes. This procedure allows us to distinguish between non-adiabatic effects coming 
from shocks and the ones from cooling. The role of the cooling, even when it could be 
important in other scenarios, is irrelevant for the simulations considered in this paper, while 
shocks play the most important role. 

In the picture describing the dynamics of the baryonic component, there are some clues 
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showing the presence of shocks. Examining the quantities sensitive to shocks, all of them 
evidence the formation of a quasi-spherical shock. This shock seems to arise around z ~ 2 
at the cluster center and moves outwards. Nevertheless, some irregular shocks could form 
at z > 2. This conclusion is supported by the behaviour of the entropy profiles (see Fig. 5), 
and the existence of shocked cells at these times (see Fig. 8). The quasi-spherical shock 
would form from the collapse of the quasi-spherical global structure , while other smaller 
shocks - with a filamentary morphology - would arise from some collapsing substructure 
and merging processes. In short, previous discussion manifest two different regimes in the 
shock formation. 

It should be noticed that the structures simulated in this paper correspond, due to the 
initial conditions, to a large Abell cluster. For this kind of clusters, gravitational collapse 
is fast and the dynamics is violent. Shocks form earlier and are stronger than in others 
smaller cluster-hke objects (< 3a). 

Some discussion on the numerical resolution of the simulations is needed. The one 
used in present paper (~ O.SMpc) is not enough to simulate the very center of the clusters 
and galaxy formation, but it suffices to study the role of the shocks in ICM. It should 
be kept in mind that HRSC techniques arc able to resolve shocks even with coarse grids. 
Nevertheless, higher resolution would be desiderable to perform more complete simulations. 
Improvements in numerical resolution will introduce smaller scales in the problem, as a 
consequence, the physics of the model should be enriched in order to describe this new 
scenario. Chemical reactions and radiative transfer should be considered. 
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A. Energy conservation properties 

As it is well known, in simulations without cooling the total energy must be conserved. 
This conservation law can be obtained after integrating on the whole computational volume 
the evolution equations, and it reads as follows (Peebles 1980): 

M±Z)^ipE + Vn=0 (Al) 

dt a 

where E is the total energy - kinetic of the gas plus kinetic of the particles plus internal for 
the gas- and W is the total gravitational energy, i.e., gas plus particles. 



The Eq. (|Al|) can be intregated respect to the scale factor giving: 

a{t2)E{t2) - a{ti)E{ti) + T^*'^ Eda = -{a{t2)W{t2)) - a{ti)W{ti)) (A2) 

Ja(ti) 

being t2 and ti two different times. The energy conservation can be tested by defining a 
quantity, R, as: 

. . ^ a(t)^(t) - a{t,)E{t,) + Eda 

-{a{t)W{t))-a{t,)W{t^)) ^ ' 



Consequently with Eq. (X2), quantity R must remain equal to unit during the evolution. 



We have tested the energy conservation for the whole code , i.e., hydrodynamical 
and particle codes, coupled through the multidimensional Poisson solver. One numerical 
simulation with the same initial condition as in the case A, and with 64'^ cells and 64'^ 
particles has been considered. Figure 12 plots the quantity R against redshift. The maxima 



-23- 



errors are around three per cent. As it is well known, the better resolution (finer grids and 
more particles), the better conservation properties. 

B. Some comments on the resolution of the hydrodynamical code 

Due to the coarse grids we have used in this paper, our results depend strongly on the 
spatial resolution. Our numerical simulations have been performed using the hydro-code 
described in Quilis et al. (1996). In that paper the authors discussed two types of 
reconstruction procedures as methods to increase the spatial resolution of the algorithm, 
i.e. MUSCL (linear) and PPM (parabolic). In order to investigate the influence of the 
numerical resolution in our galaxy cluster simulations, we have considered some numerical 
experiments. In order to do that, some simulations without cooling and with the same 
initial conditions than in A case (see Section 2) have been performed. 

The above two reconstruction procedures have been considered. Together with the 
reconstruction, we have studied the influence of the number of cells and particles used in 
the simulations. Extremely coarse grids with 64^ cells and 64^ particles, and the one used 
in the simulations in this paper, 128^ cells and 128^ particles, have been considered. 

Table 1 summarizes the main results of these experiments. We show the maxima 
of the density contrast of the gas , 5^"^, at two times. At each time, these maxima are 
presented for both grids and set of particles , i.e. 64^ and 128^, and for MUSCL and PPM 
reconstructions. As it follows from that table, with the grids considered, MUSCL does 
not work properly. PPM appears as the best procedure to increase the spatial resolution, 
in this sense, the simulations presented in this paper have been computed using PPM 
reconstruction procedure. 
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Table 1: Influence of the numerical resolution. 

z MUSCL64 PPM64 

gmax I 24.5 56.3 
gmax 11.5 96.0 



MUSCL128 PPM128 
53.2 210.5 
31.5 348.8 
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Fig. 1. — Dark matter particles distribution inside the computational boxes for six epochs for 
A simulation. Boxes are 20 comoving Mpc length per edge. Only 32^ particles are displayed 
in each box. 

Fig. 2. — Slices along z-axis of the gas density contrast for A simulation. Each slice is 20 
comoving Mpc length per edge and 0.15 comoving Mpc depth. Columns stand for fixes 
redshift. Central rows are centered at the maximum density contrast, up(low) row is located 
at -5(+5) comoving Mpc from the central slices along z-axis. At left side of each slice there 
is one palette describing the colour scale used to plot it. 

Fig. 3. — Average density contrast, for A simulation, in radial shells for gas(top) and 
DM(bottom). Six times are displayed. Radii are in comoving Mpc. 

Fig. 4. — Average radial velocities, for A simulations, in radial shells for gas(top) and 
DM(bottom). Six times are displayed. Radii are in comoving Mpc and velocities in units of 
speed of light (c). 

Fig. 5. — Average entropy in radial shells for gas in A simulations. Six times are displayed. 
Radii are in comoving Mpc. Entropy is defined as s = Zn(T/p^~^). 

Fig. 6. — From top to bottom, slices show (for the baryonic component): divergence of 
velocity in units of c(first row), entropy (second row), temperature in K (third row), and 
pressure in dyn/cm^ (fourth row). Each slice is 20 comoving Mpc length per edge and 0.15 
comoving Mpc depth. Columns stand for a fixed redshift. All shoes are centered at the 
maxima of the density contrast. At left side of each shoe there is one palette describing the 
colour scale. Results correspond to A simulation. 

Fig. 7. — Shocked cells inside the computational boxes for six times for A simulation. Boxes 
are 20 comoving Mpc length per edge. 
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Fig. 8. — Number of shocked cells (Nghc) versus redshift. Continuous(daslied) line plots 
A(B) simulation. 

Fig. 9. — Plot of logarithm of the total internal energy {Ejj) against redshift (top panel) and 
the ratio of the total internal energy (Eu) to the total kinetic (Ex) energy versus redshift 
(bottom panel). All energies are in units of ergs. Continuous (dashed-pointed) line stands 
for A(B) results, in top panel dashed line shows analytical adiabatic result. 

Fig. 10. — Plot of the logarithm of the average temperature (in K degrees) in a ball of 1 
comoving Mpc centered at the core cluster. Continuous, dashed-pointed, dashed lines stands 
for A and B simulations, and analytical adiabatic evolution, respectively. 

Fig. 11. — Plot for B simulation of the logarithm of the X-Ray luminosity (in erg/s), from 
a ball of 1 comoving Mpc, centered at the core cluster. 

Fig. 12. — Plot for the quantity R defined in Eq. (|A3|) against the redshift. Perfect energy 
conservation would imply R = 1. 
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